Synthesizing by means of mice

1 Introduction

First, we again create 1 complete dataset, which is done by means of imputing the boys dataset in mice (Van Buuren & Groothuis-Oudshoorn, 2011) with default settings, and complete the dataset. Then, based on this one true dataset, we extract the regression coefficients of a linear regression model in which weight is regressed on age and height, that are used to compare the estimates of the synthetic versions of the datasets. The true regression coefficients are equal to \(\beta_0 = -9.216\), \(\beta_{age} = 2.327\), \(\beta_{height} = 0.191\). Then, for a total of 200 iterations, we impute a matrix with the same dimensions as the boys dataset (that is, the sample size \(n = 748\) and the number of predictors \(k = 9\)). This is done once with the default settings, that is, with pmm for the continuous variables, polr for binary variables and polyreg for ordinal variables with more than two categories. Furthermore, bmi is imputed passively, since it consists of a fixed combination of height and weight bmi = (wgt / (hgt/100)^2. Furthermore, the predictor matrix is adjusted so that imputations for bmi do not flow back into predictions of hgt and wgt. Then, for every synthetic dataset, the estimates, and corresponding variance of the estimates are calculated, as well as the \(95\%~CI\) are calculated. Furthermore, an indicator is added, whether or not the confidence interval includes the true parameter estimates. Note that these confidence intervals are still based on imputations for missing data, instead of imputations for synthetic data. This is due to the fact that the calculation for the degrees of freedom for the estimates based imputing synthetic values might yield degrees of freedom that are so small that the corresponding \(95\% ~ CI\) yields boundary values of \(-\infty\) and \(\infty\).

---

2 Simulations

---

2.1 Single dataset approach

---

2.1.1 True results

broom::tidy(truemodel, conf.int=TRUE) %>%
  mutate(CIW = conf.high - conf.low) %>% 
  kable(digits = 3,
        caption = "Results of a single linear regression model on the completed dataset.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Results of a single linear regression model on the completed dataset.

term	estimate	std.error	statistic	p.value	conf.low	conf.high	CIW
(Intercept)	-9.216	2.033	-4.534	0	-13.207	-5.226	7.982
age	2.327	0.189	12.335	0	1.956	2.697	0.741
hgt	0.191	0.028	6.843	0	0.136	0.246	0.110

---

2.1.2 Use mice to synthesize and to pool

---

2.1.2.1 Pool the results of the synthetic datasets with the default mice settings.

default_fit <- syns_def %>%
  map(function(x) x %$% lm(wgt ~ age + hgt) %>% pool %>% summary) %>%
  map_dfr(function(x) {
    var      <- x[,1]
    true_est <- coefs
    est      <- x[,2]
    true_se  <- sqrt(diag(vcov(truemodel)))
    se       <- x[,3]
    df       <- x[,5]
    lower    <- est + se * qt(.025, df)
    upper    <- est + se * qt(.975, df)
    cov      <- lower < coefs & coefs < upper
  
    bind_cols(var = var, true_est = true_est, est = est, true_se = true_se, 
                               se = se, df = df, lower = lower, upper = upper, cov = cov)
    })

results_def <- default_fit %>%
  group_by(var) %>%
  summarise("True Est" = unique(true_est),
            "Syn Est"  = mean(est),
            "Bias"     = mean(est - true_est),
            "True SE"  = unique(true_se),
            "Syn SE"   = mean(se),
            "df"       = mean(df),
            "Lower"    = mean(lower),
            "Upper"    = mean(upper),
            "CIW"      = mean(upper - lower),
            "Coverage" = mean(cov))

2.1.2.2 Pool the results of the synthetic dataset with CART.

cart_fit <- syns_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt) %>% pool %>% summary) %>%
  map_dfr(function(x) {
    var      <- x[,1]
    true_est <- coefs
    est      <- x[,2]
    true_se  <- sqrt(diag(vcov(truemodel)))
    se       <- x[,3]
    df       <- x[,5]
    lower    <- est + se * qt(.025, df)
    upper    <- est + se * qt(.975, df)
    cov      <- lower < coefs & coefs < upper

    bind_cols(var = var, true_est = true_est, est = est, true_se = true_se,
              se = se, df = df, lower = lower, upper = upper, cov = cov)
    })

 results_cart <- cart_fit %>%
   group_by(var) %>%
   summarise("True Est" = unique(true_est),
             "Syn Est"  = mean(est),
             "Bias"     = mean(est - true_est),
             "True SE"  = unique(true_se),
             "Syn SE"   = mean(se),
             "df"       = mean(df),
             "Lower"    = mean(lower),
             "Upper"    = mean(upper),
             "CIW"      = mean(upper - lower),
             "Coverage" = mean(cov))

2.1.2.3 Results for both synthesization methods

bind_rows("Mice default" = results_def,
          "Mice with cart" = results_cart, .id = "Imputation method") %>%
  kable(digits = 3,
        caption = "Results of with mice default settings, and with CART settings, for the imputation of synthetic data, with the in-build mice pooling function.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Results of with mice default settings, and with CART settings, for the imputation of synthetic data, with the in-build mice pooling function.

Imputation method	var	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
Mice default	(Intercept)	-9.216	-2.160	7.057	2.033	2.941	28.084	-8.657	4.338	12.995	0.338
Mice default	age	2.327	2.948	0.621	0.189	0.285	18.624	2.304	3.591	1.287	0.490
Mice default	hgt	0.191	0.094	-0.097	0.028	0.041	22.553	0.003	0.186	0.183	0.388
Mice with cart	(Intercept)	-9.216	-7.602	1.614	2.033	2.664	50.879	-13.212	-1.992	11.220	1.000
Mice with cart	age	2.327	2.463	0.136	0.189	0.256	41.663	1.917	3.008	1.091	1.000
Mice with cart	hgt	0.191	0.169	-0.022	0.028	0.037	43.755	0.090	0.249	0.159	1.000

---

2.1.3 Use mice to synthesize with pooling rules from Drechsler

For now, we forget about the default mice imputation method, because it leads to flawed results. Therefore, we continue with the CART results only. We fit the lm model in which wgt is regressed on age and hgt.

## pool.lm.syn

## If you have a multiply imputed synthetic dataset, that is, multiple fully imputed 
## datasets that are generated by mice, using an all-ones matrix as the where matrix, 
## this function can be used to analyse the synthetic data by means of a linear 
## regression model, and to pool the results.

pool.syn <- function(mira) {
  
  if(class(mira)[1] == "mira") {
    fitlist <- mira %$% analyses
  }
  else {
    fitlist <- mira
    }
  
  m <- length(fitlist)
  
  pooled <- fitlist %>% 
    map_dfr(broom::tidy) %>%
    group_by(term) %>%
    summarise(est     = mean(estimate),
              bm      = sum((estimate - est)^2) / (m - 1),
              ubar    = mean(std.error^2),
              var_u   = (1 + 1/m) * bm - ubar,
              var     = if_else(var_u > 0, var_u, ubar),
              df      = max(1, (m - 1) * (1 - ubar / (bm + bm/m))^2),
              lower   = est - qt(.975, df) * sqrt(var),
              upper   = est + qt(.975, df) * sqrt(var), .groups = 'drop')
  pooled
}

pool2.syn <- function(mira) {
  
  if(class(mira)[1] == "mira") {
    fitlist <- mira %$% analyses
  }
  else {
    fitlist <- mira
  }
  
  m <- length(fitlist)
  
  pooled <- fitlist %>% 
    map_dfr(broom::tidy) %>%
    group_by(term) %>%
    summarise(est     = mean(estimate),
              bm      = sum((estimate - est)^2) / (m - 1),
              ubar    = mean(std.error^2),
              var_u   = (1 + 1/m) * bm - ubar,
              var     = if_else(var_u > 0, var_u, ubar),
              df      = max(m - 1, (m - 1) * (1 - ubar / (bm + bm/m))^2),
              lower   = est - qt(.975, df) * sqrt(var),
              upper   = est + qt(.975, df) * sqrt(var), .groups = 'drop')
  pooled
}

pool3.syn <- function(mira) {
  
  if(class(mira)[1] == "mira") {
    fitlist <- mira %$% analyses
  }
  else {
    fitlist <- mira
  }
  
  m <- length(fitlist)
  
  pooled <- fitlist %>% 
    map_dfr(broom::tidy) %>%
    group_by(term) %>%
    summarise(est     = mean(estimate),
              bm      = sum((estimate - est)^2) / (m - 1),
              ubar    = mean(std.error^2),
              var     = ubar + bm/m,
              df      = (m - 1) * (1 + (ubar * m)/bm),
              lower   = est - qt(.975, df) * sqrt(var),
              upper   = est + qt(.975, df) * sqrt(var), .groups = 'drop')
  pooled
}

ci_cov <- function(pooled, true_fit) {
  
  coefs <- coef(true_fit)
  vars   <- diag(vcov(true_fit))
  
  nsim <- nrow(pooled) / length(coefs)
  
  pooled %>% mutate(true_coef = rep(coefs, nsim),
                    true_var  = rep(vars, nsim),
                    cover     = lower < true_coef & true_coef < upper) %>%
    group_by(term) %>%
    summarise("True Est" = unique(true_coef),
              "Syn Est"  = mean(est),
              "Bias"     = mean(est - true_coef),
              "True SE"  = unique(sqrt(true_var)),
              "Syn SE"   = mean(sqrt(var)),
              "df"       = mean(df),
              "Lower"    = mean(lower),
              "Upper"    = mean(upper),
              "CIW"      = mean(upper - lower),
              "Coverage" = mean(cover), .groups = "drop")
}

models <- syns_cart %>% map(function(x) x %$% lm(wgt ~ age + hgt))

Then, we use the custom pooling function to pool the estimates.

pooled <- models %>% map_dfr(pool.syn)

And another custom function to summarise all results, and obtain the confidence interval coverage.

ci_cov(pooled, truemodel) %>%
  kable(digits = 3,
        caption = "Imputation results of synthesizing by means of CART, with the pooling rules by Drechsler.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Imputation results of synthesizing by means of CART, with the pooling rules by Drechsler.

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-7.602	1.614	2.033	1.811	72.945	-20.997	5.793	26.789	1
age	2.327	2.463	0.136	0.189	0.169	91.939	1.106	3.820	2.714	1
hgt	0.191	0.169	-0.022	0.028	0.025	112.809	-0.027	0.365	0.392	1

It can be seen that the coverage of the confidence interval is nearly equal to one. Since we currently only want to make inferences with regard to the sample, this is okay. Namely, if the CI coverage with regard to the sample would be \(95\%\), the CI coverage with regard to the population would be \(0.9025\%\). Inferences with regard to the population will be discussed in the next paragraph.

However, what is questionable, is the fact that the confidence interval width has increased, even though the variance has decreased. Since the confidence interval width is a function of only the standard error and the degrees of freedom, this must be due to the fact that the degrees of freedom are smaller than the degrees of freedom as returned by the mice pooling function.

mice_df <- data.frame(term = cart_fit$var, df = cart_fit$df)
syn_df <- data.frame(term = pooled$term, df = pooled$df)

df <- bind_rows("Mice" = mice_df, "Synthetic" = syn_df, .id = "Method")

ggplot(data = df, aes(x = df, color = Method, fill = Method)) +
  geom_density(alpha = .7) +
  facet_wrap(~ term, nrow = 1, scales = "free") +
  scale_color_viridis_d(begin = .30, end = .65) +
  scale_fill_viridis_d(begin = .30, end = .65) +
  theme_classic()

---

2.1.4 Different synthesizing sequence

Now we have established that the mice synthesizing algorithm works as good as the synthpop synthesizing algorithm, we can check whether changing the order in which hgt and wgt are synthesized influences the result of the algorithm.

cart_wgt_hgt %>% 
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-7.635	1.581	2.033	1.800	40.577	-20.105	4.834	24.939	0.994
age	2.327	2.459	0.132	0.189	0.165	43.528	1.171	3.747	2.576	0.996
hgt	0.191	0.170	-0.021	0.028	0.025	34.469	-0.015	0.354	0.369	0.998

---

2.2 Bootstrapped boys data

true_results <- bootstrap_boys %>%
  map(~ lm(wgt ~ age + hgt, .x)) %>%
  map_dfr(~ broom::tidy(.x, conf.int = TRUE)) %>%
  mutate(true_coef = rep(coef(truemodel), nsim),
         true_se   = rep(sqrt(diag(vcov(truemodel))), nsim),
         cover     = conf.low < true_coef & true_coef < conf.high) %T>% 
  assign("boot_true", ., envir = globalenv()) %>%
  group_by(term) %>%
  summarise("True Est" = unique(true_coef),
            "Boot Est"  = mean(estimate),
            "Bias"     = mean(estimate - true_coef),
            "True SE"  = unique(true_se),
            "Boot SE"   = mean(std.error),
            "DF"       = 745,
            "Lower"    = mean(conf.low),
            "Upper"    = mean(conf.high),
            "CIW"      = mean(conf.high - conf.low),
            "Coverage" = mean(cover)) %>%
  kable(digits = 3,
        caption = "Bootstrapped estimates of regression coefficients, with mean standard error and coverage based on the individual estimates of the standard error.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

synthetic_results <- boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %T>% 
  assign("boot_syn", ., envir = globalenv()) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3,
        caption = "Synthetic (CART) results with pooling rules by Drechsler.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

norm_results <- boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  mutate(df = NA,
         lower = est - qnorm(.975) * sqrt(var),
         upper = est + qnorm(.975) * sqrt(var)) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3, caption = "Normal CI approximation") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Bootstrapped estimates of regression coefficients, with mean standard error and coverage based on the individual estimates of the standard error.

term	True Est	Boot Est	Bias	True SE	Boot SE	DF	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-9.301	-0.084	2.033	2.025	745	-13.275	-5.326	7.949	0.934
age	2.327	2.314	-0.013	0.189	0.188	745	1.945	2.683	0.738	0.932
hgt	0.191	0.193	0.001	0.028	0.028	745	0.138	0.247	0.109	0.922

Synthetic (CART) results with pooling rules by Drechsler.

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-8.113	1.104	2.033	1.828	147.690	-19.472	3.247	22.719	0.968
age	2.327	2.418	0.091	0.189	0.167	61.161	1.289	3.547	2.258	0.964
hgt	0.191	0.176	-0.015	0.028	0.024	92.660	0.014	0.339	0.325	0.958

Normal CI approximation

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-8.113	1.104	2.033	1.828	NA	-11.696	-4.529	7.167	0.860
age	2.327	2.418	0.091	0.189	0.167	NA	2.091	2.744	0.653	0.818
hgt	0.191	0.176	-0.015	0.028	0.024	NA	0.128	0.224	0.096	0.790

We find that some bias has been introduced due to synthesizing the relationships, so we investigate this in further depth. Additionally, the coverage is somewhat too high, and the average confidence interval width has increased dramatically. Therefore, we will further investigate where this increase stems from.

plot_dat <- bind_rows("Synthetic" = boot_syn %>% 
                                    mutate(ciw = upper - lower) %>%
                                    select(term, est, ciw = ciw), 
                      "Real values" = boot_true %>% 
                                      mutate(ciw = conf.high - conf.low) %>%
                                      select(term, est = estimate, ciw = ciw),
                      .id = "Method")


ggplot(plot_dat, aes(x = est, color = Method, fill = Method)) +
  geom_density(alpha = .7) +
  facet_wrap(~ term, nrow = 1, scales = "free") +
  scale_color_viridis_d(begin = .30, end = .65) +
  scale_fill_viridis_d(begin = .30, end = .65) +
  theme_classic() +
  ggtitle("Distribution of the bootstrapped and synthetic regression coefficients.")

It can be seen that the bias is systematic, and thus not due to any extreme values that pull the average estimate to zero. Additionally, it seems that the pooled standard error indicated in the table is somewhat too low for the synthetic data. Namely, it appears in the plot that the spread in the synthetic data is roughly equal to the spread in the bootstrapped, but real, data. However, in the table above, the synthetic standard errors are much smaller than the mean over the manually calculated standard error of the bootstrapped estimates. When the standard errors are calculated over the obtained regression estimates, we find that the synthetic standard error increases much more when we calculate it over the observed synthetic estimates.

boot_se <- boot_true %>% group_by(term) %>% summarise("Mean bootstrapped SEs" = sd(estimate))
boot_syn_se <- boot_syn %>% group_by(term) %>% summarise("Mean synthetic SEs" = sd(est))

bind_cols(boot_se, boot_syn_se[,2]) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

term	Mean bootstrapped SEs	Mean synthetic SEs
(Intercept)	2.105	2.094
age	0.207	0.209
hgt	0.030	0.030

ggplot(plot_dat, aes(x = ciw, color = Method, fill = Method)) +
  geom_density(alpha = .7) +
  facet_wrap(~ term, nrow = 1, scales = "free") +
  scale_color_viridis_d(begin = .30, end = .65) +
  scale_fill_viridis_d(begin = .30, end = .65) +
  theme_classic() +
  ggtitle("Distribution of the confidence interval width (CIW) of the bootstrapped and the synthetic datasets.")

When we look at the confidence interval width, we find that the bootstrap CI is of approximately the same length in all iterations, while the CI’s of the synthetic results spread out much more.

---

2.3 Improving the bias/variance tradeoff

The fact that there is bias introduced, might be due to the fact that there is too little variance in the estimates of the synthetic values. Therefore, an approach that increases the variance of the estimates might yield estimates that display less bias. This can be achieved by decreasing the number of iterations, that is, setting the mice parameter maxit to 1, instead of 5. Additionally, the complexity parameter can be set to \(10^{-8}\) instead of \(10^{-4}\), so that more splits will be made. Also, the value of the rpart parameter minbucket, that is, the minimum number of observations in any terminal node, is set to 3 instead of 5, so that, once again, more splits can be made. By means of these changes, the within variance increases, while the between-imputations variability decreases. The results with 500 iterations on the same, singly imputed boys dataset can be seen below.

syns_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3,
        caption = "Synthesizing the same data.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Synthesizing the same data.

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-8.391	0.825	2.033	1.821	52.299	-20.543	3.761	24.305	1.000
age	2.327	2.405	0.078	0.189	0.164	29.103	1.169	3.641	2.472	0.998
hgt	0.191	0.179	-0.012	0.028	0.024	38.830	0.002	0.357	0.355	1.000

Additionally, the same approach is used to synthesize the 500 bootstrapped samples from the completed boys data.

boot_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3,
        caption = "Synthesizing bootstrap samples with adjusted imputation settings.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Synthesizing bootstrap samples with adjusted imputation settings.

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-9.080	0.137	2.033	1.952	238.879	-16.892	-1.268	15.624	0.956
age	2.327	2.332	0.006	0.189	0.178	240.207	1.484	3.180	1.696	0.952
hgt	0.191	0.190	-0.002	0.028	0.027	215.781	0.067	0.312	0.245	0.946

However, I’m not sure whether it is sensible to set the number of iterations to maxit = 1. Therefore, I additionally included an example that shows the same approach, that is, with the complexity parameter cp = 1e-08 and the parameter minbucket = 3, but with maximum number of iterations fixed at maxit = 5.

boot_cart_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3,
        caption = "Synthesizing bootstrap samples with adjusted imputation settings but with maxit = 5.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Synthesizing bootstrap samples with adjusted imputation settings but with maxit = 5.

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-9.037	0.180	2.033	1.964	250.209	-16.683	-1.390	15.293	0.956
age	2.327	2.337	0.010	0.189	0.178	218.191	1.513	3.160	1.647	0.940
hgt	0.191	0.189	-0.002	0.028	0.027	177.875	0.071	0.306	0.235	0.954

It can be seen that the bias indeed decreases, while the standard error of the estimates increases. This is due to the fact that the within-imputation variance increases, while the between-imputation variance decreases somewhat.

cart_default <- boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  group_by(term) %>%
  summarise(Within = mean(ubar), 
            Between = mean(bm), 
            Total = mean(var),
            .groups = 'drop')

cart_adj_maxit <- boot_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  group_by(term) %>%
  summarise(Within = mean(ubar), 
            Between = mean(bm), 
            Total = mean(var),
            .groups = 'drop')

cart_adj <- boot_cart_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  group_by(term) %>%
  summarise(Within = mean(ubar), 
            Between = mean(bm), 
            Total = mean(var),
            .groups = 'drop')

bind_rows("Cart default" = cart_default,
          "Cart with adjusted settings and maxit = 1" = cart_adj_maxit,
          "Cart adjusted settings" = cart_adj, .id = "Method") %>%
  kable(digits = 3,
        caption = "Within, between and total imputation variance for default CART, and adjusted CART.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Within, between and total imputation variance for default CART, and adjusted CART.

Method	term	Within	Between	Total
Cart default	(Intercept)	3.890	2.164	3.516
Cart default	age	0.033	0.022	0.030
Cart default	hgt	0.001	0.000	0.001
Cart with adjusted settings and maxit = 1	(Intercept)	4.053	1.447	3.892
Cart with adjusted settings and maxit = 1	age	0.035	0.015	0.033
Cart with adjusted settings and maxit = 1	hgt	0.001	0.000	0.001
Cart adjusted settings	(Intercept)	4.054	1.374	3.934
Cart adjusted settings	age	0.035	0.014	0.033
Cart adjusted settings	hgt	0.001	0.000	0.001

---

3 Adjusted degrees of freedom

Next to the analyses displayed previously, I adjusted the degrees of freedom to \(\text{max}(m - 1, \nu)\), as is done in Reiter and Drechsler (2010). This yields the following results for the initial analyses, that is, bootstrapping with the mice default settings, but with the method set to cart. Note that this concerns 500 bootstrap samples from the complete data.

boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool2.syn) %T>% 
  assign("boot_syn", ., envir = globalenv()) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-8.113	1.104	2.033	1.828	149.153	-12.616	-3.610	9.006	0.906
age	2.327	2.418	0.091	0.189	0.167	62.838	1.998	2.837	0.839	0.886
hgt	0.191	0.176	-0.015	0.028	0.024	94.288	0.115	0.237	0.122	0.866

Additionally, the results can be computed for the data with the smaller bias settings. That is, setting maxit = 1, cp = 1e-08 and minbucket = 3.

boot_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool2.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-9.080	0.137	2.033	1.952	239.625	-13.563	-4.597	8.966	0.940
age	2.327	2.332	0.006	0.189	0.178	241.218	1.911	2.753	0.842	0.914
hgt	0.191	0.190	-0.002	0.028	0.027	216.711	0.127	0.252	0.125	0.908

And, once again, the same analysis as the previous one, but with maxit = 5.

boot_cart_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool2.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

term	True Est	Syn Est	Bias	True SE	Syn SE	df	Lower	Upper	CIW	Coverage
(Intercept)	-9.216	-9.037	0.180	2.033	1.964	250.899	-13.523	-4.550	8.973	0.936
age	2.327	2.337	0.010	0.189	0.178	219.162	1.917	2.756	0.839	0.902
hgt	0.191	0.189	-0.002	0.028	0.027	178.770	0.127	0.251	0.125	0.930

References

Van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate imputation by chained equations in r. Journal of Statistical Software, 45(3), 1–67. Retrieved from https://www.jstatsoft.org/v45/i03/